/**
 * @file ebex_pdf.h
 *
 * @date Mar 4, 2010
 * @author seth
 *
 * @brief This file is part of EBEX State Estimator, created for the EBEX project
 *
 * This software is copyright (C) 2010 Columbia University
 *
 * EBEX State Estimator is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * EBEX State Estimator is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with EBEX State Estimator; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#ifndef EBEX_PDF_H_
#define EBEX_PDF_H_

#include <stddef.h>
#include <ebex_lookup.h>

#define _kernel_types(x,_)	\
	_(x, none)				\
	_(x, epanechnikov)		\
	_(x, uniform)			\
	_(x, triangle)			\
	_(x, gaussian)			\
	_(x, quartic)			\
	_(x, triweight)			\
	_(x, cosine)			\
	_(x, tricubic)

EBEX_LOOKUP_TABLE(kernel_type, static);

void ebex_stats_pdf(int ndata, double data[], int dstride, double weight[], int wstride, double h,
		e_kernel_type m_kernel, int m_outputnum, double *m_outputx,
		int m_outputx_stride, double *m_outputy);
void ebex_stats_cpdf(int m_numdata, const double *m_data, int m_datastride, double m_bandwidth, e_kernel_type m_kernel,
		int m_outputnum, double *m_outputx, int m_outputx_stride, double *m_outputy);
void ebex_stats_cpdfc(int m_numinputdata, double *m_data, int m_datastride,
		double m_bandwidth, e_kernel_type m_kernel, int m_numoutputdata,
		double m_outputx[], int m_outputx_stride, double *m_outputy);
double ebex_stats_bandwidth(int m_numdata, double *m_data,
		size_t m_datastride, double *m_weight,
		size_t m_weightstride, e_kernel_type m_kernel);

/* Kernel-based local regression */
void ebex_stats_smooth(
		size_t m_numinput,
		double* m_inputy, size_t m_inputy_stride,
		double *m_inputx, size_t m_inputx_stride,
		double m_bandwidth, e_kernel_type m_kernel,
		size_t m_numoutput,
		double *m_outputx, size_t m_outputx_stride,
		double *m_outputy, size_t m_outputy_stride);

double ebex_stats_mean(const size_t m_dim, const double *m_data);
double ebex_stats_dbl_mean_variance(size_t m_numdata, const double *m_data, size_t m_datastride, bool m_biased, double *m_mean);
double ebex_stats_dbl_weighted_delta_variance(size_t m_numdata, const double *m_data, const double *m_weight, bool m_biased);
double ebex_stats_median(size_t m_dim, double *m_data);
int ebex_stats_covariance(int m_num, int m_dim, const double *m_data, size_t m_major_stride, size_t m_minor_stride, bool m_biased, double *m_cov);
void ebex_stats_weighted_covariance(int m_num, int m_dim, const double *m_data, const double *m_mean, const double *m_weight,
                                   size_t m_major_stride, size_t m_minor_stride, bool m_biased, double *m_cov);
int
ebex_stats_lowess(const double *m_inputx, const double *m_inputy, size_t m_numinput,
		size_t m_windowsize,
		size_t m_numiter,
		double m_delta,
		double *m_outputy, size_t m_outputy_stride,
		double *m_robust_weight, double *m_residual);


static inline double ebex_stats_residual_epanechnikov(double m_offset, double m_bandwidth)
{
	double u = m_offset / m_bandwidth;
	return (1 - u * u);
}

static inline double ebex_stats_convolve_epanechnikov(const double *m_inputy, size_t m_inputy_stride,
                                                      const double *m_inputx, size_t m_inputx_stride,
                                                      double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth) continue;
		if (u > m_bandwidth) break;
		double val = ebex_stats_residual_epanechnikov(u, m_bandwidth);
		numerator += val * m_inputy[j * m_inputy_stride];
		denominator += val;
	}
	if (denominator > 0.0)
		retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_uniform(double m_offset __attribute__((unused)), double m_bandwidth __attribute__((unused)))
{
	return 1.0;
}

static inline double ebex_stats_convolve_uniform(const double *m_inputy, size_t m_inputy_stride,
                                                 const  double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = m_inputx[j * m_inputx_stride] - m_outputx;
		if (u < -m_bandwidth)
			continue;
		if (u > m_bandwidth)
			break;
		double val = ebex_stats_residual_uniform(u, m_bandwidth);
		numerator += val * m_inputy[j * m_inputy_stride];
		denominator += val;
	}
	if (denominator > 0.0) retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_triangle(double m_offset, double m_bandwidth)
{
	double u = m_offset / m_bandwidth;
	return fabs(u);
}

static inline double ebex_stats_convolve_triangle(const double *m_inputy, size_t m_inputy_stride,
                                                  const double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth)
			continue;
		if (u > m_bandwidth)
			break;
		double val = ebex_stats_residual_triangle(u, m_bandwidth) ;
		numerator += val * m_inputy[j * m_inputy_stride];
		denominator += val;
	}
	if (denominator > 0.0) retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_gaussian(double m_offset, double m_bandwidth)
{
	double u = m_offset / m_bandwidth;
	return exp(-0.5 * u * u);
}

static inline double ebex_stats_convolve_gaussian(const double *m_inputy, size_t m_inputy_stride,
                                                  const double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth)
			continue;
		if (u > m_bandwidth)
			break;
		double val = ebex_stats_residual_gaussian(u, m_bandwidth);
		numerator += m_inputy[j * m_inputy_stride] * val;
		denominator += val;
	}
	if (denominator > 0.0) retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_quartic(double m_offset, double m_bandwidth)
{
	double u = m_offset / m_bandwidth;
	return pow(1 - u * u, 2.0);
}

static inline double ebex_stats_convolve_quartic(const double *m_inputy, size_t m_inputy_stride,
                                                 const double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth)
			continue;
		if (u > m_bandwidth)
			break;
		double val = ebex_stats_residual_quartic(u, m_bandwidth);
		numerator += m_inputy[j * m_inputy_stride] * val;
		denominator += val;
	}
	if (denominator > 0.0) retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_triweight(double m_offset, double m_bandwidth)
{
	double u = m_offset / m_bandwidth;
	return pow(1 - u * u, 3.0);
}

static inline double ebex_stats_convolve_triweight(const double *m_inputy, size_t m_inputy_stride,
                                                   const double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth)
			continue;
		if (u > m_bandwidth)
			break;
		double val = ebex_stats_residual_triweight(u, m_bandwidth) ;
		numerator += m_inputy[j * m_inputy_stride] * val;
		denominator += val;
	}
	if (denominator > 0.0) retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_cosine(double m_offset, double m_bandwidth)
{
	return cos(M_PI_2 * m_offset / m_bandwidth);
}

static inline double ebex_stats_convolve_cosine(const double *m_inputy, size_t m_inputy_stride,
                                                const double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth)
			continue;
		if (u > m_bandwidth)
			break;
		double val = ebex_stats_residual_cosine(u, m_bandwidth) ;
		numerator += m_inputy[j * m_inputy_stride] * val;
		denominator += val;
	}
	if (denominator > 0.0) retval = numerator / denominator;

	return retval;
}

static inline double ebex_stats_residual_tricubic(double m_offset, double m_bandwidth)
{
	return pow(1 - pow(fabs(m_offset / m_bandwidth),3.0), 3.0);
}

static inline double ebex_stats_convolve_tricubic(const double *m_inputy, size_t m_inputy_stride,
                                                  const double *m_inputx, size_t m_inputx_stride,
                                                    double m_bandwidth, double m_outputx, size_t m_numinput)
{
	double retval = 0.0;
	double numerator = 0.0;
	double denominator = 0.0;
	for (size_t j = 0; j < m_numinput; j++)
	{
		double u = (m_inputx[j * m_inputx_stride] - m_outputx);
		if (u < -m_bandwidth) continue;
		if (u > m_bandwidth) break;
		double val = ebex_stats_residual_tricubic(u, m_bandwidth) ;
		numerator += val * m_inputy[j * m_inputy_stride];
		denominator += val;
	}
	if (denominator > 0.0)
		retval = numerator / denominator;

	return retval;
}
#endif /* EBEX_PDF_H_ */
